TK!2
=v
#1
:n

:s
b
:c
>>> Numerical Differentiation <<<
#2
:n
fun
:v
'lnā
:s
i
:c
name of function to differentiate
#3
:n
h
:v
.00001
:s
i
:c
increment
#4
:n
x
:v
1.25
:s
i
:c
independent variable
#5
:n
y3
:v
-.2274535334093053
:s
o
:c
1st derivative, 3-point formula
#6
:n
y5
:v
-.2274535333670936
:s
o
:c
1st derivative, 5-point formula
#7
:n
dy%dx
:v
1.197310028899778
:s
o
:c
2nd derivative, 3-point formula
#8
:n
dy%dxr
:v
1.197337402836229
:s
o
:c
2nd derivative, recursive formula
=u
=l
=f
#1
:n
diff13
:c
First derivative, 3-point formula
:t
Procedure
:i
fun,x,h
:o
z
#2
:n
diff15
:c
First derivative, 5-point formula
:t
Procedure
:i
fun,x,h
:o
z
#3
:n
diff23
:c
Second derivative, 3-point formula
:t
Procedure
:i
fun,x,h
:o
z
#4
:n
rdiff
:c
Recursive Differentiation Procedure
:t
Procedure
:p
fun,h
:i
n,x
:o
z
#5
:n
lnā
:c
ln(gamma(x)) for evaluating psi function
:t
Rule
:a
x
:v
z
#1
:n
diff13
:z
" Notation: fun  - name of the differentiated function
"           x    - value of the independent variable
"           h    - increment

z:= (apply(fun,x+h)-apply(fun,x-h))/(2*h)
/E
E
#2
:n
diff15
:z
" Notation: fun  - name of the differentiated function
"           x    - value of the independent variable
"           h    - increment

z:=(apply(fun,x-2*h)-8*(apply(fun,x-h)-apply(fun,x+h))-apply(fun,x+2*h))/h/12
/E
E
#3
:n
diff23
:z
" Notation: fun  - name of the differentiated function
"           x    - value of the independent variable
"           h    - increment

z:= (apply(fun,x+h)-2*apply(fun,x)+apply(fun,x-h))/h^2
/E
E
#4
:n
rdiff
:z
if n>1 then z:= (rdiff(n-1,x+h)-rdiff(n-1,x-h))/(2*h)
if n=1 then z:= (apply(fun,x+h)-apply(fun,x-h))/(2*h)

" This function returns an approximate value of the nth derivative of fun(x)
" at  x ;  the parameter variable  h  carries the increment in the "three-
" point" numerical differentiation formula used here.
" The function f(x) must be defined by a rule, procedure or built-in
" function the name of which is carried from the Variable Sheet by another
" parameter variable  fun .
" The round-off and truncation errors in numerical differentiation are
" generally hard to control.  As a rule of thumb, h=1e-6 is a good choice for
" the 1st derivative when the order of magnitude of f(x) is 1; for the nth
" derivative one should use h=10^(7-n).  Larger  h  will increase the trun-
" cation error, but reduce the round-off error that becomes critical with
" higher order derivatives.
/E
E
#5
:n
lnā
:z
z = ln(gamma(x))

" Note: dz/dx = dln(gamma(x))/dx = psi(x)  is so-called psi or digamma f-n
"       d2z/dx2 = dpsi(x)/dx  is trigamma function
/E
E
=r
#1
:r
y3 = diff13(fun,x,h)       " This model exemplifies the usage of procedures
#2
:r
y5 = diff15(fun,x,h)       " based on several numerical differentiation for-
#3
:r
dy%dx = diff23(fun,x,h)    " mulas (see the Function Sheet for the summary
#4
:r
dy%dxr = rdiff(2,x)        " and the Function Subsheets for explanation).
#5
:r
" The above solution relates to the so-called psi or digamma and trigamma
:s
C
#6
:r
" functions (see f-n lnā). The results compare well with -.2274535334 and
:s
C
#7
:r
" 1.1973291545 from NBS Handbook of Mathematical Functions, 1972, p.267.
:s
C
#8
:r
" Another experiment: assign fun='sin and x=pi()/4 and compare the results
:s
C
#9
:r
" with sqrt(.5)=.7071067811865476 .
:s
C
=p
=t
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
Yes
:a
Yes
:r
Yes
:h
No
:t
1
:n
Yes
:d
LPT1
%Tv,f,0,5,0,1,0
%Or,f,0,1,0,1,0
